//////////////////////////////////////////////////////////////////////////
// boundtopology - class to hold information about topology of boundary
// By Antonio Miranda - Fev 2013
//////////////////////////////////////////////////////////////////////////

#include <math.h>

#include "boundtopology.h"
#include "../topology/tetratopology.h"
#include "../ctopology/topology.hpp"


#define EDGE_COMMON   -1
#define EDGE_BOUND     1
#define EDGE_CRACK     2
#define EDGE_MARK      3

#define SURF_COMMON   -1
#define SURF_BOUND     1
#define SURF_CRACK     2
#define SURF_MARK      3

#define NODE_COMMON    0
#define NODE_CRACK     2
#define NODE_MARK      3

static double angleEdge (int i, cTopology *topol);
static void createMeshFromSeedElement( int id, cTopology *m_currMesh, vector<int> &elemVec );
static void getNodesCrack ( int idnode, int curr_edge, cTopology *m_currMesh, vector<int> &ndVec);

//
//////////////////////////////////////////////////////////////////////////
boundTopology::boundTopology( tetraTopology *tetratopol )
{
  m_tetraTopol = NULL;

  // basic checks
  if (tetratopol == NULL)
    return;

  m_tetraTopol = tetratopol;
  m_surfTopol  = NULL;

  // build boundary surface topology
  buildSurfToppology();

  // identify surfaces
  identifySurfaces();

  identifyCracks ();
}

//
//////////////////////////////////////////////////////////////////////////
boundTopology::~boundTopology(void)
{
  delete m_surfTopol;
}


//
//////////////////////////////////////////////////////////////////////////
void boundTopology::getCrackNodes( vector<int> &nodes )
{
  vector<int> temp;

  for (size_t i = 0; i < m_crackFront.size(); ++i)
  {
    for (size_t j = 0; j < m_crackFront[i].size(); ++j)
       nodes.push_back(m_invMapNode[m_crackFront[i][j]]);
  }

}


//
//////////////////////////////////////////////////////////////////////////
void boundTopology::buildSurfToppology()
{
  int i, j, id_node;

  // get max node id and reserve vector space to mapping
  int max_id = m_tetraTopol->maxIdNode();
  m_mapNode.resize(max_id+1);
  for (i = 0; i < max_id+1; ++i)
    m_mapNode[i] = -1;

  // get nodes and boundary faces
  vector <node *> nodes = m_tetraTopol->getNodes ();
  vector <triface *> bound = m_tetraTopol->getBound ();

  // variables to insert to cTopology
  double *coords; 
  int    nnodes = 0, *conn;
  int    nfaces = (int) bound.size();
  conn = new int[nfaces*4];

  // build element connectivity
  for (i = 0; i < nfaces; ++i)
  {
    conn[i*4+0] = 3; // triangle element
    for (j = 0; j < 3; ++j)
    {
      id_node = bound[i]->vert[j];
      if (m_mapNode[id_node] == -1)
      {
        m_mapNode[id_node] = nnodes;
        ++nnodes;
      }
      conn[i*4+j+1] = m_mapNode[id_node];
    }
  }

  // inverse mapping node
  m_invMapNode.resize(nnodes);

  // nodes
  coords = new double[nnodes*3];
  for (i = 0; i < max_id+1; ++i)
  {
    if (m_mapNode[i] != -1)
    {
      for (j = 0; j < 3; ++j)
        coords[m_mapNode[i]*3+j] = nodes[i]->coord[j];

      m_invMapNode[m_mapNode[i]] = i;
    }
  }

  // create surface topology
  if (m_surfTopol == NULL)
    m_surfTopol = new cTopology;
  else
    m_surfTopol->Clear();

  // InsertMesh
  m_surfTopol->InsertMesh(nnodes, nfaces, coords, conn, NULL);

//   m_surfTopol->writeMatlab ("boundary.m", "wt");

  delete []conn;
  delete []coords;

}

//
//////////////////////////////////////////////////////////////////////////
void boundTopology::identifySurfaces()
{
  int i, j, nadj, adj1;
  double angle;

  // clear all crack edges 
  m_crackEdges.clear();

  // set all edges flags
  int nedges = m_surfTopol->NumEdges();
  for (i = 0; i < nedges; ++i)
  {
    m_surfTopol->SetEdgeInfo(i, EDGE_COMMON);
    angle = angleEdge(i, m_surfTopol);

    if (angle == 360.0)
    {
      m_crackEdges.push_back (i);
      m_surfTopol->SetEdgeInfo(i, EDGE_CRACK);
    }
    else if (angle <= 135.0)
      m_surfTopol->SetEdgeInfo(i, EDGE_BOUND);
  }

  // set all faces as SURF_COMMON
  int n_faces = m_surfTopol->NumElems();
  for (i = 0; i < n_faces; ++i)
    m_surfTopol->SetElemInfo(i, SURF_COMMON);

  // identify surfaces crack
  int n_crackEdges = (int) m_crackEdges.size();
  for (i = 0; i < n_crackEdges; ++i)
  {
    // for all faces adjacent to crack edge 
    nadj = m_surfTopol->NumAdjElemToEdge (m_crackEdges[i]);
    for (j = 0; j < nadj; ++j)
    {
      vector<int> elemVec;
      adj1 = m_surfTopol->AdjElemToEdge (m_crackEdges[i], j);

      // from a face (element), find other faces
      createMeshFromSeedElement (adj1, m_surfTopol, elemVec);

      // hold list of faces
      if (elemVec.size() > 0)
      {
        markSurfaces(elemVec, SURF_CRACK);
        m_crackFaces.push_back(elemVec);
      }
    }
  }

}

//
//////////////////////////////////////////////////////////////////////////
void boundTopology::identifyCracks()
{
  int ei, ej;
  int n_edges = (int) m_crackEdges.size();
  if (n_edges == 0)
    return;

  // set all nodes on crack tip as NODE_CRACK
  for (int i = 0; i < n_edges; ++i)
  {
    m_surfTopol->GetEdge (m_crackEdges[i], &ei, &ej);
    m_surfTopol->SetNodeInfo (ei, NODE_CRACK);
    m_surfTopol->SetNodeInfo (ej, NODE_CRACK);
  }

  // set all nodes on crack tip as NODE_CRACK
  for (int i = 0; i < n_edges; ++i)
  {
    m_surfTopol->GetEdge (m_crackEdges[i], &ei, &ej);

    // get nodes sequencely from i and j
    vector<int> nodesVeci, nodesVecj;
    vector<matrix4x4> tmp;

    getNodesCrack(ei, m_crackEdges[i], m_surfTopol, nodesVeci);
    getNodesCrack(ej, m_crackEdges[i], m_surfTopol, nodesVecj);

    // join nodes
    if (nodesVeci.size() > 0 || nodesVecj.size() > 0)
    {
      // add nodes to the end of nodeVeci
      for (size_t j = 0; j < nodesVecj.size(); ++j)
        nodesVeci.push_back(nodesVecj[j]);

      // hold front
      m_crackFront.push_back(nodesVeci);

      // reserve transform matrix
      tmp.resize(nodesVeci.size());
      m_crackTransf.push_back(tmp);
    }
  }

  // compute matrix transformation for all crack points
  // makeCrackTransf ();

}

// 
//////////////////////////////////////////////////////////////////////////
void boundTopology::markSurfaces (vector<int> ids, int flag)
{
    size_t n = ids.size();
    for (size_t i = 0; i < n; ++i)
    {
        m_surfTopol->SetElemInfo(ids[i], flag);
    }
}


// 
//////////////////////////////////////////////////////////////////////////
void boundTopology::makeCrackTransf ( )
{
  int pti, ptj;
  
  // number of cracks
  size_t n_cracks = m_crackFront.size();
  
  // for all cracks
  for (size_t i = 0; i < n_cracks; ++i) 
  {
    // number of points on crack tip
    size_t n_pts = m_crackFront[i].size();
    
    // check for internal crack
    bool internal_crack = false;
    int  first_node = m_crackFront[i][0];
    int  end_node   = m_crackFront[i][n_pts-1];
    int  edge       = m_surfTopol->GetEdge (first_node, end_node);
    if (edge != -1)
      internal_crack = true;
    
    // for all points on crack tip
    for (size_t j = 0; j < n_pts; ++j) 
    {
      
      // initial point j = 0
      if (j == 0)
      {
        if (internal_crack) 
        {
          pti = end_node;
          ptj = m_crackFront[i][1];
        }
        else
        {
          pti = m_crackFront[i][0];
          ptj = m_crackFront[i][1];
        }
      }
      else if (j == n_pts-1)  // end point 
      {
        if (internal_crack) 
        {
          pti = m_crackFront[i][n_pts-2];
          ptj = first_node;
        }
        else
        {
          pti = m_crackFront[i][n_pts-2];
          ptj = m_crackFront[i][n_pts-1];
        }
      }
      else // intermediary nodes 
      {
        pti = m_crackFront[i][j-1];
        ptj = m_crackFront[i][j+1];
      }
      
      // compute crack transformation for this position
      computeCrackTransf(i, j, pti, ptj);    

    }
  }
  
}

// 
//////////////////////////////////////////////////////////////////////////
void boundTopology::computeCrackTransf (int id_crack, int id_pos, int id_tang0, int id_tang1)
{
  Point3D pt0, pt1, z_axis, pt_center, curr_pt, pt_tmp;
  double  tmp[3], center[3];
  
  // current node
  int curr_node = m_crackFront[id_crack][id_pos];
  m_surfTopol->GetCoordNode(curr_node, tmp);
  curr_pt.setXYZ(tmp);
  
  // get first node to compute tangent
  m_surfTopol->GetCoordNode(id_tang0, tmp);
  pt0.setXYZ(tmp);

  // get second node to compute tangent
  m_surfTopol->GetCoordNode(id_tang1, tmp);
  pt1.setXYZ(tmp);
  
  // tangent
  z_axis = pt1 - pt0;

  
  // adjacent nodes to crack node
  int ne, *elem, n_faces = 0;
  m_surfTopol->AdjElemToNode(curr_node, &ne, &elem);
  
  // get a point on crack surface
  pt_center.setXYZ(0.0, 0.0, 0.0);
  for (int i = 0; i < ne; ++i) 
  {
    if (m_surfTopol->GetElemInfo(elem[i]) != SURF_CRACK)
      continue;
    
    // get center
    m_surfTopol->GetElemCenter(elem[i], center);
    pt_center.setX(pt_center.getX()+center[0]);
    pt_center.setY(pt_center.getY()+center[1]);
    pt_center.setZ(pt_center.getZ()+center[2]);
    
    ++n_faces;
  }
  pt_center /= n_faces;
  delete []elem;
  
  // vector from curr_node to median position of faces
  pt_tmp = pt_center - curr_pt;
  
  // axis
  Point3D x_axis, y_axis;
  y_axis = pt_tmp.crossProd(z_axis);
  x_axis = y_axis.crossProd(z_axis);
  x_axis.normalize();
  y_axis.normalize();
  z_axis.normalize();
  
  // set transformation matrix
  m_crackTransf[id_crack][id_pos].set(x_axis.getX(), y_axis.getX(), z_axis.getX(), curr_pt.getX(),
                                      x_axis.getY(), y_axis.getY(), z_axis.getY(), curr_pt.getY(),
                                      x_axis.getZ(), y_axis.getZ(), z_axis.getZ(), curr_pt.getZ(),
                                      0.0, 0.0, 0.0, 1.0);
  
  //Point3D test(1.0, 1.0, 0.0), test2;
  //test2 = m_crackTransf[id_crack][id_pos].map(test);
  
  
}


// compute angle between adjacent faces
//////////////////////////////////////////////////////////////////////////
double angleEdge( int edge, cTopology *topol )
{
  int i, nadj, adj1, adj2;
  double cross[3], dot, v1[3], v2[3], len;
  double len1, len2;

  static double pi = fabs (acos (-1.0));

  nadj = topol->NumAdjElemToEdge (edge);
  adj1 = topol->AdjElemToEdge (edge, 0);
  adj2 = topol->AdjElemToEdge (edge, 1);

  // normal
  topol->GetElemNorm (adj1, v1);
  topol->GetElemNorm (adj2, v2);
  len1 = sqrt (v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
  len2 = sqrt (v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]);

  if (len1 == 0.0 || len1 == 0.0)
    return 180.0;

  /* normalize vector */
  for (i = 0; i < 3; i++)
  {
    v1[i] /= len1;
    v2[i] /= len1;
  }

  dot   = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
  cross[0] = v1[1]*v2[2]-v1[2]*v2[1];
  cross[1] = v1[2]*v2[0]-v1[0]*v2[2];
  cross[2] = v1[0]*v2[1]-v1[1]*v2[0];
  len = sqrt (cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2]);

  // deal with 360 and 180
  if (len == 0.0)
  {
    if (dot >= 0.0) 
      return 180.0;
    else 
      return 360.0;
  }

  if (dot > 1.0)
    dot = 1.0;
  if (dot < -1.0)
    dot = -1.0;

  return (acos (dot) * 180.0 / pi);
}


//
//////////////////////////////////////////////////////////////////////////
void createMeshFromSeedElement( int id, cTopology *m_currMesh, vector<int> &elemVec )
{
  int nedges, i, j;
  int numNodes = m_currMesh->NumNodes();
  int numElements = m_currMesh->NumElems();

  list <int> edgesList;

  // add first edges to list
  nedges = m_currMesh->NumAdjEdgeToElem (id);
  for (i = 0; i < nedges; ++i)
  {
    int id_edge = m_currMesh->AdjEdgeToElem (id, i);
    edgesList.push_back(id_edge);
  }

  // for all edges until find the end
  list <int> ::iterator itr;
  for (itr = edgesList.begin(); itr != edgesList.end(); itr++)
  {
    int idEdge = *itr;

    // check for edge
    int flagEdge = m_currMesh->GetEdgeInfo(idEdge);
    if (flagEdge != EDGE_COMMON)
      continue;        // go to next edge

    // set edge flag as 1
    m_currMesh->SetEdgeInfo(idEdge, EDGE_MARK);

    // look for adjacent elements to edges
    int numElem = m_currMesh->NumAdjElemToEdge (idEdge);
    for (i = 0; i < numElem; ++i)
    {
      int idElem = m_currMesh->AdjElemToEdge (idEdge, i);
      int flagElem = m_currMesh->GetElemInfo(idElem);
      if (flagElem != SURF_COMMON)
        continue;   // go to next element

      // set element flag as 1
      m_currMesh->SetElemInfo(idElem, SURF_MARK);
      elemVec.push_back (idElem);

      // add adjacent edges to element
      int nAdjedges = m_currMesh->NumAdjEdgeToElem (idElem);
      for (j = 0; j < nAdjedges; ++j)
      {
        int id_edge = m_currMesh->AdjEdgeToElem (idElem, j);
        int flag_edge = m_currMesh->GetEdgeInfo(id_edge); 
        if (flag_edge == EDGE_COMMON)
          edgesList.push_back(id_edge);
      }
    }
  }

}

//
//////////////////////////////////////////////////////////////////////////
void getNodesCrack( int idnode, int curr_edge, cTopology *topol, vector<int> &ndVec )
{
  int ei, ej;

  if (topol->GetNodeInfo(idnode) == NODE_MARK)
    return;

  int curr_node = idnode;
  int next_node;

  do 
  {
    // add current node
    ndVec.push_back(curr_node);
    topol->SetNodeInfo(curr_node, NODE_MARK);

    // get next edge different to curr_edge
    int n_edges = topol->NumAdjEdgeToNode (curr_node);
    next_node = -1;
    for (int i = 0; i < n_edges; ++i)
    {
      int id_edge = topol->AdjEdgeToNode (curr_node, i);

      // only edge_crack and different of curr_edge
      if (id_edge == curr_edge || topol->GetEdgeInfo(id_edge) != EDGE_CRACK)
        continue;

      topol->GetEdge (id_edge, &ei, &ej);
      if (ei == curr_node)
        next_node = ej;
      else
        next_node = ei;

      if (topol->GetNodeInfo(next_node) == NODE_MARK)
      {
        next_node = -1;
        continue;
      }
      else
      {
        curr_node = next_node;
        break;
      }
    }

  } while (next_node != -1);


}
